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Abstract. 

This  paper  reports  an  investigation  into  the  use  of  CLD  methods  for  modelling  the  combined  fluid 
flow  and  heat  transfer  in  a  generic  combustor  heatshield  geometry  typical  of  current  aeroengine 
designs.  The  study  describes  the  development  and  testing  of  a  conjugate  heat  transfer  methodology,  to 
allow  the  CLD  code  to  predict  the  heat  exchange  in  both  fluid  and  solid  regions  of  the  solution  domain 
simultaneously;  special  attention  is  given  to  continuity  of  heat  flux  at  the  fluid/solid  interface.  The 
heatshield  cooling  air  is  supplied  via  rows  of  jets  that  flow  through  holes  in  the  backplate  and  impinge 
on  the  heatshield  backside.  An  interesting  phenomenon  associated  with  these  jets  is  investigated.  It  is 
noticed  that  under  some  conditions  of  heatshield  flow,  jet  impingement  leads  to  the  model  predicting  a 
minimum  in  heat  transfer  coefficient,  rather  than  a  maximum,  at  the  impingement  point.  It  is  shown, 
by  predicting  the  measured  Nusselt  number  distribution  in  a  single  impinging  jet  experiment,  that  this 
is  a  correct  simulation,  associated  with  the  combined  effects  of  low  impingement  height  and  low 
turbulence  levels  in  the  jet  cores.  The  presence  of  pedestals  on  the  impingement  surface  removes  this 
effect  by  creating  extra  turbulence  and  enhancing  conductive  paths  for  heat  transfer.  A  model  for 
including  the  pedestals  in  both  flow  and  heat  exchange  processes  as  a  sub-grid-scale  model  in  the  CLD 
simulation  is  described.  An  illustrative  calculation  of  the  performance  of  the  overall  model  for  a 
realistic  heatshield  geometry  is  provided,  indicating  the  predicted  spatial  valuation  of  temperature  on 
heatshield  surfaces. 

1.  Introduction. 

The  estimation  of  the  temperature  field  throughout  the  solid  material  of  a  gas-turbine  combustor  is  an 
important  element  of  the  overall  design  process.  Combustor  wall  temperatures  are  an  essential  factor 
which  exert  a  strong  influence  on  the  induced  thermal  stress,  material  properties,  and  fatigue/creep 
calculations  that  limit  the  life  of  the  combustor.  The  choice,  and  detailed  design  parameters,  of  the 
cooling  air  strategy  are  driven  to  a  large  extent  by  desired  maximum  metal  surface  temperatures,  to  be 
achieved  using  the  minimum  amount  of  cooling  air.  These  considerations  control,  for  example,  the 
number  and  location  of  small  film  cooling  slots,  size  of  effusion  patches  containing  thousands  of  small 
holes,  use  of  cooling  tiles  with  multiple  pedestal  arrays  etc. 

Accurate  and  early  knowledge  of  the  temperature  field  is  particularly  important  in  the  case  of 
combustor  components  which  are  highly  stressed  regions,  or  are  subject  to  high  heat  transfer  rates, 
such  as  the  combustor  heatshield.  The  heatshield  sits  in  the  head  of  the  combustor,  surrounding  the 
fuel  injector,  and  acts  as  a  thermal  shield  against  the  extremely  hot,  radiating  gases  in  the  combustor 
primary  zone.  It  must  be  cooled  by  air  drawn  through  the  cowl  of  the  combustor,  and  minimising  this 
cooling  air,  subject  to  allowable  metal  temperatures  and  heatshield  stress  fields,  is  a  difficult  design 
task.  This  task  is  further  complicated  by  the  complex  geometry  of  this  part  of  the  combustor,  see  Lig. 
1.  The  cooling  airflows  into,  through  (i.e.  within  the  backplate  and  heatshield  ‘sandwich’),  and  out  of 
the  heatshield  front  face,  follow  complex  3D  flow  paths.  To  maximise  heat  transfer  rates,  many 
heatshield  designs  make  use  of  heat  transfer  augmentation  devices  such  as  hu  ge  numbers  of  pin-fin 
pedestal  arrays  introduced  between  the  heatshield  and  its  backplate.  It  is  therefore  perhaps  not 
surprising  that  the  current  design  of  this  component  is,  to  a  large  extent,  empirically  driven.  The 
development  of  an  accurate  numerical  model  for  prediction  of  the  heat  transfer  processes  within 
typical  heatshield  designs  would  be  an  extremely  valuable  tool,  and  this  is  therefore  the  focus  of  the 
work  described  in  the  present  paper. 
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(a) 


Figure  1.  (a)  combustor  geometry,  (b)  heatshield  cooling  flow,  (c)  heatshield  geometry. 


It  is  extremely  rare  that  combustor  CFD  calculations  solve  simultaneously  for  flow,  reaction  and 
conjugate  heat  transfer  (i.e.  convective/radiative  heat  input  from  the  hot  gases,  convective  cooling  and 
conductive  heat  transfer  through  solid  components).  The  usual  approach  (e.g.  Krebs  et  al  [1])  is  to  use 
measured,  or  estimated,  wall  temperatures  as  boundary  conditions  for  the  CFD  predictions.  Baldi  et  al 
[2]  have  provided  perhaps  the  most  complete  description  of  a  combined  thermo-fluids  and 
mechanical/materials  stress  analysis  approach  to  life  prediction  of  combustor  liners,  using  a  hybrid 
CFD/FEM  model.  The  combustor  geometry  considered  was,  however,  unrepresentative  in  its 
simplicity,  and  more  emphasis  was  placed  on  the  mechanical  modelling  of  liner  creep  than  on  the 
accuracy  of  the  heat  transfer. 

Crocker  et  al  [3]  present  CFD  calculations  of  a  complete  combustor  geometry  from  compressor  exit  to 
turbine  inlet  (including  flow  predictions  both  external  and  internal  to  the  combustor  liner).  These 
authors  did  include  both  convective  and  conduction  heat  transfer  and  hence  produced  liner  wall 
temperature  predictions.  The  combustor  components  did  not,  however,  include  any  of  the  small-scale 
features  mentioned  above  as  this  would  have  lead  to  prohibitively  large  meshes  if  these  features  were 
to  be  resolved  by  the  grid.  Full  coupling  of  numerical  fluid  flow  and  conjugate  heat  transfer  models  is 
more  common  in  the  field  of  turbine  blade  cooling,  e.g.  Bohn  et  al  [4].  Flere  similar  challenges  occur 
as  described  above  in  the  combustor  context,  i.e.  complex  geometry,  small  scale  heat  exchange 
features  (pedestals  are  often  used  in  the  thin  trailing  edge  regions  of  blades),  and  minimising  of 
cooling  flows.  Even  in  this  application,  however,  internal  pin-fin  arrays  are  often  excluded  on  grid 
resolution  grounds,  even  though  these  contribute  significantly  to  the  heat  transfer.  Recently,  Luff  and 
McGuirk  [5]  have  shown  how  a  sub-grid-scale  modelling  approach  may  be  adopted  to  take  account  of 
the  pedestal  heat  transfer  without  unnecessarily  enlarging  the  computational  task. 
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The  present  paper  continues  reporting  of  the  work  begun  in  Ref.  [5].  The  same  numerical  techniques 
are  used  for  heatshield  and  pedestal  modelling.  The  opportunity  is  taken  here  to  demonstrate  the 
effectiveness  of  the  approach  adopted  for  conjugate  heat  transfer.  In  Ref.  [5]  a  curious  observation  was 
made  when  the  model  was  applied  to  a  full  heatshield  calculation.  This  concerned  the  appearance  of 
local  minima  in  heat  transfer  coefficients  at  cooling  air  jet  impingement  points,  most  noticeable  in  the 
absence  of  pedestals.  The  prediction  of  this  phenomenon  is  investigated  and  validated  here  before  the 
model  is  again  applied  to  the  heatshield  geometry  with  a  small  improvement  made  in  the  treatment  of 
the  effusion-type  cooling  holes  which  form  a  cooling  film  over  the  heatshield  front  surface. 


The  following  section  provides  details  of  the  mathematical  model  and  numerical  scheme  used.  This  is 
followed  by  an  investigation  of  jet  impingement  heat  transfer  for  the  short  impingement  height  condi¬ 
tions  relevant  to  heatshield  flows.  The  3D  coupled  fluid  flow/conjugate  heat  transfer  model  is  than  ap¬ 
plied  to  a  combustor  heatshield  geometry  generically  similar  to  that  found  in  current  Rolls-Royce 
engine  designs.  Since  the  emphasis  in  the  present  paper  is  on  description  of  the  computational 
methodology,  only  ‘isolated’  heatshield  calculations  are  presented,  with  simple  but  representative 
boundary  conditions  selected  to  represent  combustor  primary  zone  heat  input. 


2.  Mathematical  Model  and  Numerical  Techniques. 

2.1  Equations  governing  fluid  flow. 

The  continuity  and  momentum  equations  are  written  in  Cartesian  tensor  form.  Steady  mean  flow  is  as¬ 
sumed.  The  approximation  is  also  made  that  there  is  no  internal  heat  generation  and  the  temperature 
changes  are  sufficiently  small  that  the  fluid  properties  (in  particular  the  density)  may  be  assumed  con¬ 
stant.  (NB  it  is  very  easy  to  relax  this  assumption  and  allow  the  fluid  density  to  be  calculated  from  an 
ideal  gas  relation,  but  the  simplification  is  made  here  as  illustrative  calculations  only  are  presented). 
Under  these  conditions,  the  temperature  becomes  a  passive  scalar,  but  it  is  important  to  emphasise  that 
the  energy  equation  is  eventually  to  be  solved  in  both  fluid  and  solid  regions  of  the  chosen  solution 
domain. 
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Upper  case  letters  indicate  time-averaged  mean  values  and  lower  case  letters  fluctuations  about  the 
mean.  U,  is  the  velocity,  P  the  pressure  and  T  the  temperature;  p  is  the  fluid  density,  cp  is  specific  heat 
and  vmoi  and  Xmoi  represent  the  diffusive  properties  for  momentum  (kinematic  viscosity)  and  heat  trans¬ 
fer  (conductivity)  respectively.  The  standard  high  Reynolds  number  two-equation  k-£  model  and  asso¬ 
ciated  wall  functions  (Launder  and  Spalding  [6])  are  used;  the  Reynolds  stresses  and  turbulent  heat 
fluxes  in  Eqns.  (2)  and  (3)  are  thus  modelled  via: 


UiUj  =Vt 


3U;  3Uj 

- L  + - - 

dx:  dx ■ 


—  k5y 

3  J 


(4) 


—  X,  9T 


■ujt  = 


pcp  dxj 


(5) 


vt=Cu  — 


(6) 


Ref.  [6]  provides  details  on  the  full  transport  equations  for  k  and  £. 
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If  the  energy  equation  is  considered  as  a  typical  example,  insertion  of  the  turbulence  closure  scheme 
leads  to: 


^eff  3T 
pep  dxj 


(7) 


^eff  ^mol  + 
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Where  the  molecular-  and  turbulent  Prandtl  numbers  (<rT  and  at,i)  are  taken  as  0.7  and  0.9  respectively. 


In  order  to  treat  complex-shaped  domains,  the  Cartesian  versions  of  all  transport  equations  are  trans¬ 
formed  into  a  general,  curvilinear,  non-orthogonal  (q,  q,  Q  co-ordinate  system,  defined  by  a  grid 
generation  process  to  fit  the  boundary  shape  of  the  desired  geometry.  The  energy  equation  (7),  for 
example,  then  assumes  the  form: 
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Where  U,  V  and  W  now  represent  the  mass-flux  carrying  velocity  components  and  qij  are  the  metric 
components  of  the  Jacobian  tensor  J  describing  the  transformation  (for  details  see  Ref.  [7]). 


The  transformed  equations  are  solved  using  a  multi-block,  structured  grid,  finite-volume  formulation 
using  an  implicit  pressure-correction  method  (for  details  see  McGuirk  and  Spencer  [11]).  Boundary- 
fitted  non-orthogonal  grids  were  constructed  using  a  grid-generator  written  to  allow  a  CAD  solid 
model  interface  (Eccles  [12]).  The  solid  model  of  the  generic  heatshield  geometry  shown  in  Fig.  1(c) 
was  used  to  create  a  mesh  in  both  fluid  and  solid  regions  of  the  heatshield;  examples  of  the  mesh, 
which  consisted  of  3  blocks  and  around  700,000  cells  are  given  in  Ref.  [5]. 

2.2  Conjugate  heat  transfer. 

The  energy  equation  (9)  is  to  be  solved  in  both  solid  and  fluid  regions  to  allow  conjugate  heat  transfer. 
In  the  solid  region  the  velocity  components  are  zero,  so  only  the  right  hand  side  of  Eqn.  (9)  remains,  X, 
is  clearly  zero  and  Xmoi  is  now  the  solid  material  conductivity  Xs.  It  has  been  emphasised  by  other 
authors  (e.g.  Refs.  [4]  and  [8])  that  there  are  distinct  advantages  to  being  able  to  use  identical 
numerical  treatments  for  solving  the  energy  equation  in  linked  fluid/solid  regions,  but  that  this  then 
requires  careful  attention  to  be  paid  to  the  heat  flux  evaluation  at  the  fluid/solid  interface.  The  fact  that 
the  diffusion  coefficient  changes  sharply  across  a  fluid/solid  interface  (sometimes  by  orders  of 
magnitude)  means  that  discontinuous  temperature  gradients  result.  Conventional  discretisation 
practices  for  cell-centred  finite-volume  schemes  are  no  longer  accurate  (see  Patankar  [9],  who 
recommends  using  harmonic  rather  than  arithmetic  averaging  for  cell  face  conductivity).  On  highly 
non-orthogonal  grids,  as  typically  result  from  the  geometrical  complexity  of  heatshields,  the 
formulation  of  the  cross-derivative  terms  in  Eqn.  (9)  also  needs  careful  treatment.  A  numerical  scheme 
has  been  developed  that  successfully  addresses  these  issues.  This  is  described  in  detail  in  Luff  [10], 
but  an  indication  of  the  concepts  used,  and  evidence  of  their  effectiveness,  is  provided  here. 
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Figure  2.  Heat  transfer  between  two  neighbouring  finite-volume  cells. 

The  methodology  adopted  is  first  described  for  two  neighbouring  solid  regions  with  different  thermal 
conductivities,  and  then  extended  to  allow  fluid/solid  interfaces.  This  is  done  to  enable  the  concept 
behind  the  method  to  be  better  understood,  and  because  it  facilitated  the  validation  of  the  method  on  a 
test  problem  for  which  both  an  analytical  solution  and  other  authors’  numerical  results  were  available. 
The  method  is  described  using  a  2D  example,  but  has  been  generalised  and  applied  to  the  fully  3D 
case. 

Consider  the  diagram  in  Fig.  2;  this  illustrates  a  typical  scenario  for  two  neighbouring  finite-volume 
cells  (centred  at  P  and  N  respectively)  for  a  non-orthogonal  grid.  It  is  assumed  that  the  material  in  cells 
P  and  N  can  have  very  different  values  of  thermal  conductivity  (it  is  similarly  assumed  below  that 
fluid/solid  interfaces  are  arranged  in  the  grid-generation  step  to  coincide  with  cell  boundaries).  The 
heat  flow  rate  across  the  north  face  of  cell  P  may  be  written: 

qnp  =^P—  (Tb  -Tp.)  (10) 

ynP 

where  a  and  ynP  are  geometric  quantities  calculated  from  grid  information.  The  temperature  at  point  P’ 
may  be  calculated  from  (again  dP  and  1P  are  grid  calculated  quantities): 

Tp'=^(Te-Tw)+TP  (11) 

ip 

This  gives  rise  to  the  following  expression  for  the  heat  flow  rate  at  the  north  face  of  cell  P,  and  a 
similar  expression  for  the  heat  flow  rate  at  the  south  face  of  cell  N  may  be  derived: 

q„P=^p— (tb-tp)-3p-^Mt,-t„)  d2) 

ynP  y^p 

q.N  =K  — (Tb  -Tn  )+  K  -^7-(T„  -Tnw )  (13) 


Since  energy  must  be  conserved  at  cell  faces,  even  for  discontinuous  temperature  gradients, 
qnP  =  —  qsN  and  an  equation  for  the  interface  temperature  TB  may  be  obtained: 
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This  is  now  used  to  eliminate  TB  from  the  cell  face  heat  flow  rate  expressions  to  give,  for  example: 
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Note  that  this  expression  still  contains  unknown  temperatures  at  cell  faces  (e.g.  Te,  Tn,  Tne  etc.).  These 
are  normally  obtained  by  linear  interpolation,  but  in  regions  where  steep  and  discontinuous 
temperature  gradients  occur  this  can  introduce  significant  errors.  By  examination  of  Eqn.  (14)  for  the 
boundary  temperature,  it  can  be  seen  that  this  is  a  weighted  average  of  the  two  cell  temperatures  either 
side  of  the  boundary,  plus  some  extra  terms  to  account  for  differences  in  grid  skewness  on  either  side 
of  the  shared  boundary.  Neglecting  these  extra  terms  (since  the  cell  face  is  fitted  to  the  boundary 
shape,  grid  skewness  is  unlikely  to  change  dramatically  on  either  side  of  the  boundary),  a  similarly 
weighted-average  formulation  can  be  used  to  provide  relations  for  calculating  the  unknown  face 
temperatures,  for  example: 
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The  presence  of  the  cell  thermal  conductivities  in  these  relations  means  that  a  good  approximation  is 
obtained  despite  the  discontinuous  temperature  gradients  which  can  exist  if  the  cell  face  corresponds 
to  a  material  conductivity  change  or  a  solid/fluid  interface.  The  same  method  may  be  used  throughout 
the  solution  domain,  since  if  conductivity  is  constant,  it  reduces  to  linear  interpolation. 


In  order  to  validate  this  method,  the  test  case  for  conjugate  heat  transfer  recommended  by  Kelkar  [8] 
has  been  calculated.  This  corresponds  to  heat  conduction  in  a  composite  slab  with  internal  heat 
generation;  the  problem  possesses  an  analytical  solution  and  also  allows  non-orthogonality  errors  to  be 
assessed.  Details  of  the  geometry  and  the  computational  domain  (with  solution  domain  grid  lines 
schematically  indicated)  are  shown  in  Fig.  3. 


An  analytical  solution  for  the  non-dimensional  temperature  field  0(X,Y)  for  the  case  of  a  specified, 
spatially  varying  volumetric  heat  source  is  provided  in  Ref.  [8]  for  a  fixed  conductivity  ratio  between 
the  two  regions  shown,  C  =  A,2/A,i.  A  numerical  solution  is  found  by  solving  the  heat-conduction-only 
version  of  the  energy  equation  over  a  computational  domain  formed  by  a  parallelogram  embedded 
symmetrically  across  the  junction  of  the  composite  slab;  the  shape  of  the  computational  domain  is 
selected  to  force  a  non-orthogonal  grid  with  significant  grid  skewness.  Solutions  are  obtained  for  a 
range  of  grids  (5x10,  10x20,  20x40)  with  the  boundary  temperatures  on  the  computational  domain  and 
the  heat  source  at  internal  grid  nodes  fixed  from  the  exact  analytical  solution.  The  normalised  average 
error  between  the  computed  solution  0Cy  at  any  grid  node  (i,j)  and  the  exact  solution  0,,  is  calculated 
for  each  grid: 


E  = 


e?j-6y 


E0ij 


(17) 


Where  the  summation  is  taken  over  all  internal  grid  nodes.  Fig.  4  shows  the  variation  of  the 
normalised  error  for  the  3  grids  and  for  various  values  of  conductivity  ratio  C. 
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Several  comments  can  be  made  on  these  results.  Firstly,  for  C  =  1  the  gradient  of  the  decrease  in  error 
shows  that  the  scheme  is  second-order  accurate.  As  C  increases,  the  errors  grow  and  the  formal 
second-  order  accuracy  is  lost  due  to  the  neglect  of  the  extra  terms  in  the  interpolation.  For  very  high  C 
values,  second-order  accuracy  returns  since  the  importance  of  these  extra  terms  diminishes.  These 
observations  are  analogous  to  those  made  by  Kelkar  et  al  [8]  in  their  prediction  of  this  problem;  it 
should  be  noted,  however,  that  the  present  results  for  the  error  levels  shown  in  Fig.  4  are  about  50%  of 
those  given  in  Ref.  [8]  for  the  same  grid  density.  This  is  due  to  the  improved  interpolation  scheme  for 
the  cross-derivative  terms  suggested  above  in  Eqn.  (16). 


1.0 - -i - 1.0 


Figure  3.  Geometry  and  solution 

domain  for  conjugate  heat 
transfer  test  case. 
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Figure  4.  Normalised  error  in  conjugate 
heat  transfer  test  case  for 
various  grids  and  C  values. 


If  convective  fluid  flow  is  present  in  one  control  volume  and  solid  material  conduction  in  the  other 
(e.g.  at  the  P/N  cell  interface  in  Fig.2),  the  method  is  easily  adapted.  For  the  high  Reynolds  number 
turbulence  closure  treatment  adopted  here,  a  log-law-based  wall-function  approach  is  used  to  specify 
one  of  the  heat  flows;  for  example,  the  heat  flow  at  the  north  face  of  cell  P  is  calculated  via: 


pcpkp5C 


0.5  0.25 


LnP 


In 


''Eypk“C^ 
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V 


mol 


+p. 


(Tb-TpF-Mt.-t.) 

A  D 


(18) 


Flere  kP  represents  the  turbulence  energy  at  the  near  wall  grid  node,  and  CM,  K,  E  and  Pf  are  constants 
taken  from  the  standard  high  Reynolds  number  turbulence  model  (Ref.  [6]).  This  is  equivalent  to  the 
expression  used  in  the  conduction  only  case  (Eqn.  (12)),  but  with  a  modified  or  ‘effective’  value  used 
for  the  ‘cell’  conductivity  A,pf ,  viz.: 


\  eff 

A  p 


ynPpcpk«-5C“-25 


'  t,T 


-In 


^Eypk“crJ  ' 


mol 


+  Pf 


the  remainder  of  the  analysis  thus  proceeds  as  before  but  using  this  effective  value  for  A,P. 


(19) 


2.3  Pedestal  sub-grid-scale  model. 

To  augment  heat  transfer  between  the  hot  front  surface  of  the  heatshield  and  the  cooling  air  flowing  in 
between  the  heatshield  and  backplate,  use  is  commonly  made  of  pin-fin  type  pedestals.  The  large 
number  but  small  size  of  these  would  make  it  prohibitively  difficult  to  resolve  the  flow  around  each 
pin  and  the  heat  transfer  from  each  pin  accurately  in  the  finite-volume  model.  Luff  and  McGuirk  [5] 
have  shown  how  the  main  effects  of  the  pedestals  (in  terms  of  extra  pressure  drop,  extra  turbulence 
generation  and  extra  heat  transfer)  may  be  accounted  for  adequately  using  a  sub-grid-scale  model.  This 
consists  of  adding  extra  terms  into  the  momentum  and  turbulence  equations  and  solving  an  extra 
energy  equation  to  determine  the  local  pedestal  temperature  (purely  conduction,  but  with  anisotropic 
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conductivity,  allowing  conduction  only  in  the  direction  parallel  to  the  pedestal  axis).  For  example,  the 
pedestal  pressure  drop,  in  regions  of  the  mesh  which  are  identified  to  the  code  as  containing  pedestals, 
is  introduced  on  a  cell-  by-cell  basis  as  an  extra  pressure  gradient  in  the  momentum  equations: 


R  = 


3P 

.  dx.  . 

\  1  /extra 


=  0.5  —  pf 

s 


ped 


v.. 


V 

,max  max 


(20) 


Flere  Si  is  a  geometric  parameter  of  the  pedestal  geometry  (streamwise  pitch),  Vmax  is  the  component  of 
cell  velocity  normal  to  the  pedestals  (taking  account  of  blockage),  and  fpeci  is  an  empirical  friction 
factor  which  describes  the  extra  pressure  loss  per  row  of  pedestals.  This  latter  factor  is  taken  from 
experimental  work  on  flow  over  tube  banks,  where  pressure  loss  data  has  been  previously  correlated 
by  Zukauskas  [13]  in  terms  of  the  Reynolds  number  of  flow  over  the  pedestals,  i.e.: 


r  _ 

1  ped 


Reped 


(21) 


X,  a,  and  b  are  empirical  constants  dependant  on  pedestal  array  geometry  and  layout.  Similarly,  an 
extra  source  term  is  included  in  the  fluid  energy  equation  (and  an  equivalent  term  of  opposite  sign  in 
the  pedestal  energy  equation)  to  account  for  pedestal-to-fluid  heat  transfer,  and  this  is  again 
empirically  described  using  a  Nusselt  number  correlation,  also  taken  from  Zukauskas  [13]: 
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(22) 


This  sub-grid-scale  model  was  validated  in  Ref.  [5]  against  experimental  pressure  loss  and  heat 
transfer  data  taken  from  Metzger  et  al  [14].  The  approach  was  shown  to  perform  very  well,  in  spite  of 
the  fact  that  the  Zukauskas  [13]  empirical  correlations  were  taken  for  long  tubes  rather  than  short  pins, 
and  from  experiments  that  were  far  removed  from  the  geometrical  scenarios  typical  of  heatshield 
pedestals. 


3.  Impinging  jet  heat  transfer. 

The  cooling  airflow  in  the  heatshield  application  is  metered  through  a  series  of  round  holes  in  the 
backplate,  creating  rows  of  jets  impinging  on  the  back  surface  of  the  heatshield  (see  Fig.  1).  The 
diameter  of  these  holes,  and  the  backplate/heatshield  separation,  are  such  that  these  jets  are  at  very 
small  impingement  heights,  Fl/D.  In  earlier  work  on  combustor  heatshield  heat  transfer,  Luff  and 
McGuirk  [5]  noted  the  anomalous  result  that,  in  the  absence  of  pedestals,  local  minima  in  heat  transfer 
coefficients  were  predicted  beneath  these  impinging  jets,  rather  than  the  expected  maxima.  This  aspect 
of  heatshield  heat  transfer  is  further  examined  here. 


Impinging  jet  heat  transfer  has  been  much  studied  in  the  past,  and  has  been  noted  as  a  difficult 
problem  for  eddy  viscosity  turbulence  models  (see,  for  example,  Craft  et  al  [15]).  These  authors 
reported  that,  for  the  case  of  a  single  round  impinging  jet,  the  linear  eddy  viscosity  k  -  £  turbulence 
model  leads  to  spuriously  high  rates  of  turbulence  energy  generation  near  the  impingement  point,  and 
consequently  very  high  heat  transfer  rates,  with  Nusselt  numbers  around  twice  those  measured 
experimentally  (Ref.  [16]).  They  recommended  use  of  a  Reynolds  Stress  Transport  model  to  resolve 
this  problem.  Durbin  et  al  [17]  have  also  studied  the  same  case  and  developed  a  three-equation 
turbulence  model  to  achieve  satisfactory  predictions.  The  same  author  previously  (Ref.  [18]) 
developed  a  modification  to  the  basic  linear  k  -  £  model  (a  time  scale  limiter)  to  remove  the 
overproduction  of  k  in  stagnation  flows.  Luff  and  McGuirk  [5]  proposed  that  the  low  Nusselt  numbers 
at  impingement  in  their  heatshield  calculations  were  related  to  previously  reported  experimentally 
observed  minima  in  impingement  cooling  flows  at  low  values  of  Fl/D  and  with  low  turbulence  levels 
in  the  jet  core  (see  Ref.  [19]).  To  test  the  validity  of  this  proposal,  for  the  particular  turbulence  model 
used  here  (and  in  Ref.  [5]),  the  case  of  a  single  round  impinging  jet,  as  studied  computationally  in 
Refs.  [15]  and  [17]  and  experimentally  in  Ref.  [16],  is  predicted  here. 
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Fig.  5  presents  the  predicted  streamline  pattern  for  the  case  selected  for  study,  which  represents  the  Re 
=  23,000,  H/D  =  2.0  case  taken  from  the  measurements  reported  in  Ref.  [16].  Following  the  practice  of 
other  authors  (Refs.  [15],  [17])  the  calculations  are  started  upstream  of  the  jet  exit  to  allow  some  up¬ 
stream  distortion  of  the  jet  conditions  by  the  pressure  field  at  impingement.  Other  boundary  conditions 
are  also  as  used  by  other  authors’,  i.e.  fully-developed  velocity  and  turbulence  profiles  at  entry  to  the 
jet  pipe,  and  constant  static  pressure  boundaries  on  the  top  and  right  hand  sides  of  the  solution  domain 
(representing  essentially  entrainment  and  outflow  boundaries  as  seen  in  Fig.  5). 


Figure  5.  Streamline  pattern,  single  impinging  Figure  6.  Turbulence  energy  contours, 
jet  single  impinging  jet 

The  impingement  process  leads  to  the  formation  of  a  strong  wall-jet  flow,  whose  thickness  grows 
away  from  impingement;  the  large  k  levels  created  by  the  high  strain  rates  near  stagnation  are  also 
very  plain  to  see  (Fig.  6),  with  the  predicted  maximum  k  value  being  7.5%  of  the  square  of  the  mean 
jet  velocity.  Note  the  position  of  a  secondary  maximum  in  k  just  away  from  impingement.  It  is  known 
that  the  heat  transfer  problem  is  computationally  much  more  demanding  than  flowfield  prediction,  so  a 
series  of  grid  refinement  tests  have  been  carried  out.  Fig.  7  presents  the  predicted  Nusselt  number 
distribution  along  the  impingement  surface  for  two  grids,  where  both  axial  and  radial  mesh  refinement 
has  been  conducted.  It  is  clear  that  the  present  solutions  are  certainly  numerically  accurate  for  heat 
transfer  parameters. 

It  is  very  significant  in  Fig.  7  that  not  only  are  the  present  predictions  numerically  accurate,  but  they 
seem  to  represent  better  agreement  with  the  experimental  data  than  has  previously  been  reported. 
Although  the  stagnation  point  Nusselt  number  is  overpredicted,  the  100%  error  obtained  in  Refs.  [15] 
and  [17]  is  not  observed  here.  The  reason  for  this  is  that  the  present  results  have  been  obtained  with 
the  high-Re  wall  function  approach,  rather  than  attempting  to  use  a  low-Re  turbulence  model  and 
predict  the  flow  all  the  way  down  to  the  wall  as  was  the  case  in  Refs.  [15]  and  [17].  The  y+  value  at  the 
first  cell  next  to  the  wall  was  kept  greater  than  30  for  all  the  grids  shown  in  Fig. 7.  It  is  clear  that  most 
of  the  overprediction  of  Nu  seen  in  earlier  work  is  a  result  of  the  low  Re  model  reacting  adversely  to 
the  enhanced  k  production  near  stagnation.  When  high-Re-model-based  wall  functions  are  used, 
although  k  is  still  overpredicted,  this  does  not  lead  to  a  large  error  in  Nu.  The  peak  Nu  value  is  actually 
slightly  off-axis,  and  there  is  also  the  appearance  of  a  further  local  increase  in  Nu  in  the  region  of  the 
secondary  maximum  in  k  mentioned  above.  The  use  of  the  time-scale  limiter  recommended  by  Durbin 
[18]  is  shown  in  Fig.  8.  This  has  the  desired  effect  of  reducing  the  peak  levels  of  k  at  stagnation 
(maximum  value  now  only  5.2%  of  jet  velocity  squared),  and  this  brings  down  the  peak  Nu  levels,  to 
values  actually  slightly  below  those  measured.  The  effect  of  the  limiter  is  only  active  in  the  first  jet 
diameter  along  the  impingement  wall. 

With  the  basic  model  now  validated  as  producing  reasonable  results  for  the  experimental  case 
considered  in  Ref.  [16],  a  computational  experiment  was  now  performed,  with  the  resulting  Nusselt 
number  distribution  also  included  in  Fig.  8.  A  change  is  made  in  the  jet  nozzle  conditions;  rather  than 
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using  fully-developed  profiles,  a  flat  velocity  profile  and  uniform,  low  turbulence  entry  conditions 
were  used  (0.5%  turbulence  intensity  and  a  length  scale  of  D/10).  It  is  very  clear  to  see  how,  at  this 
low  impingement  height,  the  impingement  point  Nusselt  number  now  displays  a  local  minimum 
(reduction  by  a  factor  of  3  to  a  level  only  50%  of  the  off-axis  peak  value),  as  observed  in  previous 
experiments  (Ref.  [19])  and  full  heatshield  predictions  at  similar  values  of  H/D  (Ref.  [5]). 
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Figure  7.  Nusselt  number  predictions 
(grid  refinement). 


Figure  8.  Nusselt  number  predictions 
(time  scale  limiter  and  jet 
condition  effects). 


4.  Full  heatshield  heat  transfer. 

The  complete  model,  as  described  in  Section  2.  above,  has  been  applied  to  a  generic  heatshield  heat 
transfer  problem.  The  geometry  used  is  described  fully  in  Ref.  [5]  and  is  shown  in  Fig.  1;  symmetry 
has  been  assumed  about  the  heatshield  centre  plane.  The  cooling  air  enters  through  discrete  feed  holes 
in  the  backplate,  with  separate  rows  feeding  the  centre  section  and  upper  and  lower  sections.  The  air 
exits  from  the  latter  two  sections,  after  flowing  over  pedestal  air  ays,  as  inner  and  outer  cooling  ring 
flows  at  the  top  and  bottom  of  the  main  combustor.  The  airflow  in  the  centre  section  exits  mainly 
through  the  heatshield  front  face  via  3  annular  rings  of  small  holes.  In  Ref.  [5]  these  were  represented 
for  simplicity  as  3  circular  slots  with  the  same  geometric  area  as  the  multiple  cooling  holes.  In  the 
present  calculations  the  opportunity  has  been  taken  to  improve  the  modelling  of  this  aspect.  If  the  exit 
flows  are  treated  as  circular  slots,  it  is  believed  that  a  reasonable  representation  is  achieved  of  the  flow 
conditions;  it  is  obvious,  however,  that  the  use  of  circular  slots  provides  a  barrier  to  radial  heat 
conduction  paths,  which  are  clearly  present  in  the  real  multiple  hole  geometry.  In  the  present 
calculations,  a  similar  approach  to  that  used  to  model  pedestals  at  the  sub-grid-scale  level  has  been 
adopted  to  represent  the  individual  metal  regions  between  holes  as  ‘pin-fins’  in  the  slot  regions,  to 
allow  radial  heat  conduction  via  the  sub-grid-scale  source  terms  containing  the  pin-fin  temperature, 
calculated  via  an  extra  energy  equation. 

Numerical  treatment  and  boundary  conditions  have  otherwise  reproduced  those  used  in  Ref.  [5];  a  uni¬ 
form  total  pressure  at  all  inlets,  and  fixed  static  pressure  at  all  exits  establishes  the  mass  flow  through 
the  heatshield.  Fleat  transfer  boundary  conditions  correspond  to  a  fixed  temperature  on  the  back 
surface  of  the  backplate  (and  for  all  inlet  airflows),  at  a  level  Tin  typical  of  compressor  exit  flow.  The 
thermal  boundary  condition  on  the  front  face  of  the  heatshield  was  fixed  as  a  uniform  heat  flux  of  1 
MW/rn2;  it  would  be  relatively  easy  to  improve  this  condition  by  specifying  a  heat  flux  consistent  with 
a  convective  and  radiant  heat  input  obtained  from  a  separate  combustor  calculation,  but  it  is  believed 
that  the  simplified  condition  chosen  is  representative  and  sufficient  to  illustrate  the  methodology 
described  here. 

Fig  9  indicates  the  flow  patterns  obtained  from  the  present  model.  The  impinging  and  interacting  jets 
in  both  the  centre  and  inner/outer  portions  of  the  heatshield  are  a  very  distinctive  feature,  with  the  low 
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H/  D  value  now  being  very  obvious.  The  effect  of  the  presence  of  the  pedestals  on  the  flow  pattern  is 
mainly  to  smooth  out  the  regions  of  recirculation  between  backplate  and  heatshield. 

Figure  10  shows  the  predicted  temperature  field  on  the  hotside  of  the  heatshield.  Fig.  10(left) 
represents  predictions  ignoring  the  radial  conduction  paths  in  the  exit  slot  region  and  Fig.  lO(right) 
shows  the  effect  of  modelling  these  in  the  manner  described  above. 


Figure  9.  Velocity  vectors  in  heatshield 
flow  showing  impingement 


Figure  10.  Heatshield  hotside  temperature  contours 
(left)  without  radial  conduction  through 
slots  and  (right)  allowing  for  radial 
conduction 


It  is  noticeable  that  the  high  temperature  regions  correspond  to  regions  with  no  pedestals  or  where 
local  increases  in  metal  thickness  near  mounting  bosses  and  structural  ribs  create  greater  resistance  to 
heat  transfer.  The  allowance  for  radial  conduction  through  the  exit  slot  regions  has  not  changed  the 
overall  temperature  pattern,  and  it  has  reduced  the  peak  temperatures  by  only  around  1%. 

5.  Conclusions 

The  development  and  testing  of  a  computational  model  for  combustor  heatshield  heat  transfer  has  been 
described.  Due  to  the  complex  geometry  of  typical  heatshields,  it  is  inevitable  that  the  3D  grids  used  to 
calculate  the  cooling  flow  through  the  component  will  contain  regions  of  significant  non¬ 
orthogonality.  This  makes  the  calculation  of  conjugate  heat  transfer  through  both  fluid  and  solid 
regions  difficult,  unless  special  consideration  is  given  to  the  solid/fluid  interface  treatment.  A  method 
hat  can  cope  accurately  with  such  grid  distributions  was  described,  and  evidence  provided  of  its 
suitability  by  calculating  conduction  through  a  composite  material  with  highly  skewed  grids  and  large 
property  variations.  Heatshield  flows  contain  multiple  impinging  jets  of  cooling  air  and  this  aspect  of 
the  problem  was  also  investigated  in  some  detail.  For  a  single  impinging  jet,  the  use  of  the  standard 
high  Re-based  k-£  turbulence  model  and  log-law-based  wall-functions  was  shown  to  provide  better 
agreement  with  stagnation  point  Nusselt  number  distributions  than  was  previously  thought  possible. 
The  conditions  appropriate  to  heatshield  flows  of  low  impingement  heights  and  low  jet  core  turbulence 
levels  were  shown  to  lead  to  minima  in  heat  transfer  coefficients  at  impingement.  Finally,  illustrations 
were  provided  of  the  application  of  the  complete  model,  including  a  sub-grid-scale  model  for 
pedestals,  to  a  realistic  combustor  heatshield  geometry.  Temperature  distributions  on  the  hotside  were 
shown  to  provide  detailed  information  about  the  location  and  level  of  hot  spots.  The  model  described 
is  believed  therefore  to  have  significant  value  for  the  combustor  heatshield  design  process  and  is  a 
useful  first  step  towards  the  more  accurate  modelling  of  combustor  component  stress  fields  and 
eventually  combustor  life  prediction. 
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Question: 

How  you  determined  heat  transfer  from  the  combustor  hot  side  to  the  heat  shield? 

Answer: 

Constant  heat  input  of  1  MW/rn  2  is  prescribed  on  the  hot  side  of  the  heatshield.  This  is  a 
simplification  but  for  methodology  development  it  is  sufficient.  It  will  be  possible  in  future  to  either 
more  accurately  prescribe  heat  input,  or  to  extend  the  calculation  domain  into  the  combustor,  so  that 
convective  heat  transfer  is  predicted  on  the  hot  side  whilst  radiative  heat  input  is  prescribed. 
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